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Abstract 

Molecular  dynamics  simulations  are  used  to  demonstrate  for  the  first  time  that 
systems  about  4  nm  thick  containing  1200  to  1600  water  molecules  and  NaCl 
at  IM  to  3M  concentrations  exhibit  the  main  components  of  electric  double 
layers  at  charged  metal  surfaces.  In  particular  for  the  system  chosen  there  are 
regions  clearly  identifiable  as  a  bulk  electrolyte  zone,  a  diffuse  layer  that  screens 
the  charge  on  the  electrode  and  a  layer  of  oriented  water  localized  next  to  the 
electrode.  The  width  of  the  diffuse  layer  increases  with  decrease  in  salt  con¬ 
centration  from  0.2  nm  at  3  M  to  O.S  nm  at  1  M.  All  the  calculations  are  based 
on  the  simple  point  charge  model  (SPCE)  water  model  and  the  immersed 
electrode  approximation. 
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INTRODUCTTON 


Electric  double  layers  are  some  of  the  basic  ofgan.^  jopted  in  electrochemical  and  bi¬ 

ological  systems  to  shield  Helds  arising  from  layers  and  an  's  of  charge  in  contiguous  structures. 
It  is  therefore  quite  important  to  understand  their  properties  usinfT  simple  models  without  intro¬ 
ducing  further  approximations.  In  this  letter  we  describe  the  first  irolecular  dynamics  simulation 
of  a  complete  electric  double  layer  and  detail  it's  structure  from  the  charged  surface  of  the  metal 
electrode  all  the  way  to  the  bulk  electrolyte  zone,  in  terms  of  water  .uid  ion  time  independent 
probability  distribution  functions  averaged  parallel  to  the  metal  surface.  These  calculations 
provide  a  consistent  microscopic  picture  of  ions  and  water  throughout  a  dcuble  layer  including 
the  species  next  to  the  charged  surface,  in  the  'diffuse  layer'  and  in  the  bulk  zone.  The  effect 
of  finite  size  of  ions  and  water  are  clearly  evident,  as  is  the  effect  of  the  large  electric  field  oi 
the  orientation  of  surface  water  molecules. 

In  the  calculations  reported  here  the  simulation  cell  has  edge  length  3.724  nm,  with  a  flat  metal 
plate  at  z  =  1.862  nm  and  a  flat  restraining  wail  at  z  =  -1.862  nm.  The  cell  contains  between 
1400  to  1600  water  molecules  and  the  NaQ  concentration  is  IM,  2M  or  3M.  The  electric  charge 
density  on  the  flat  metal  plate  is  either  -9  or  -18  Cm*^. 

The  calculations  show  that  a  diffuse  layer  is  formed  between  the  metal  and  the  region  of  the  cell 
that  acts  like  a  bulk  electrolyte  zone.  The  ions  in  the  diffuse  layer  shield  the  metal  charge  so 
that  there  is  no  field  in  the  ^Ik  electrolyte  zone  which  is  locally  approximately  neutral.  In  all 
previous  papers  in  which  molecular  dynamics  or  Monte  Carlo  simulations  of  the  adsorption  of 
ions  on  electrodes  were  considered,  there  was  no  demonstration  of  the  simultaneous  existence 
of  a  diffuse  layer  and  a  bulk  region.  In  odier  words  the  simulations  were  not  on  a  scale  to  permit 
adsorption  in  the  presence  of  an  established  double  layer.  We  believe  this  is  a  prerequisite  for 
all  calculations  to  correctly  describe  the  local  and  global  properties  of  the  electric  double  layer. 
In  all  the  calculations  reported  here  we  use  the  SPCE  water  model  *  (three  charged  (qj^  = 
0.4238e)  mass  points,  bond  angle  109.5°,  OH  bond  loigth  0.1  nm,  inside  a  Lermard-Jmies  sphere 
with  radius  O  =  0.317  nm  and  well  depth  e=  0.650  kJ/mole)  and  associated  parameters  for 
NaCl2 

In  an  attempt  to  reduce  the  amount  of  computer  time  we  use  the  immersed  electrode  modeP'^ 
for  a  film  of  electrolyte  on  a  charged  metal  surface.  This  approach  is  useful  because  it  reduces 
the  number  of  water  molecules  in  the  calculation,  and  because  there  is  only  one  metal  surface 
(always  shown  on  the  rhs  in  the  figures)  there  are  no  multiple  images  in  the  calculation  of 
electric  fields  and  potentials.  In  the  immersed  electrode  model  the  charge  on  the  metal  is  deflned 
by  and  completely  screened  by  the 

excess  charge  in  the  aqueous  subphase.  It  therefore  mimics  what  happens  in  larger  systems. 
In  all  previous  simulations  the  number  of  ions  was  too  small  for  there  to  be  an  identifiable  bulk 
region  This  turns  out  not  to  be  the  case  in  the  present  set  of  calculations. 

The  metal  was  modelled  by  two  linearly  superimposed  potentials.  Pauli  repulsion  and  dispersive 
attractive  interactions  were  nwxlelled  by  a  9-3  potential,  of  the  type  used  by  Lee  et  al*.  The 
atom-suif»:e  interaction  parameters  describing  interactiem  with  nonconduction  electrons  were 
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chosen  to  be  the  same  as  those  used  by  Lee  at  al^.  The  interactions  between  ions  and  water  and 
the  conduction  electrons  was  modelled  by  their  respective  image  potentials. 

In  the  calculations  described  here  the  image  plane  and  origin  plane  of  the  9-3  potential  were 
coincident.  This  was  tantamount  to  choosing  the  image  plane  and  the  nuclear  plane  of  the  metal 
surface  to  be  coincident.  This  was  acceptable  in  our  scheme  because  the  Lennard-Jcaies  core 
parameters  o  are  all  large  and  the  'thickness'  of  the  repulsive  wall  is  also  large  (ca.  0.247  nm). 
The  calculation  of  the  image  fields  was  done  very  accurately  using  a  code  based  on  the  fast 
multipole  method  (fmm)  and  using  the  physically  correct  boundary  conditions  for  a  system  with 
one  metal  surface.  Neither  the  direct  long  range  coulomb  interaction  nor  the  image  interaction 
were  cut  off  at  a  finite  distances.  The  fast  multipole  method  developed  by  Greengard  and 
Rokhlin  is  of  order  N  in  the  number  of  charged  particles  in  the  system.  It  is  the  only  viable 
method  for  large  simulations  like  the  one  described  here. 

Finally  some  comments  on  screening  of  charged  surfaces.  At  high  salt  concentrations  the  region 
with  excess  ionic  charge  is  microscopically  small.  A  rough  estimate  of  the  thickness  of  the  zone 
in  which  screening  occurs,  valid  for  dilute  solutions  (<  0. 1  M),  is  given  by  =  ic '  the  inverse 
of  the  Debye-HUckel  screening  constant*^ 

Here  e  is  the  macroscopic  dielectric  constant  of  the  solvent  (approximately  80  for  water),  v  the 
valence  of  the  ion  (always  magnitude  one  in  this  paper),  and  n^  is  the  bulk  concentration  of  the 
ions.  Typical  values  of  d  are:  9.6  nm  for  0.001  M,  3.1  nm  for  0.01  M,  0.96  nm  for  0.1  M,  0.31 
nm  for  1.0  M,  and  0.18  nm  for  3  M  NaCl  solutions.  Obviously  in  high  salt  concentrations  the 
screening  should  be  more  efficient  and  the  screening  length  smaller,  but  since  the  numbers  are 
smaller  than  the  diameter  of  a  water  molecule  it  is  not  clear  if  they  have  any  real  meaning  out 
of  the  plane.  At  high  concentrations  there  are  in  fact  a  myriad  of  difficulties  with  simple 
Gouy-Chapman  theory and  many  modifications  to  the  theory  have  been  proposed*^.  There 
are  three  main  problems:  first  the  dielectric  constant  of  water,  second  the  finite  size  of  ions  and 
solvent  which  represent  the  lower  length  scale,  and  third  correlated  motion  amongst  ions  and 
water.  For  example  there  is  no  reason  to  believe  that  the  appropriate  value  of  e,  the  dielectric 
constant  of  water,  is  80  near  a  surface  or  in  a  high  electric  field  where  dielectric  saturation  can 
occur.  This  subject  has  been  discussed  many  times  in  the  electrochemical  literature^'^.  Therefore 
we  can  only  take  the  theoretical  values  of  screening  lengths  for  concentrated  solutions  as  an 
rough  measure  of  double  layer  thickness,  and  should  be  especially  cautious  when  d  is  smaller 
than  the  size  of  an  ion  or  solvent  molecule. 
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SCREENING  IN  SPCE  NaCI  ELECTROLYTE 

This  section  describes  the  interaction  of  hydrated  Na*^  cations  and  Cl'  anions  in  solvent  composed 
of  SPCE  water  molecules  with  a  charged  metal  surface.  Each  simulation  was  run  for  one 
nanosecond  or  longer,  and  the  instantaneous  spatial  positions  of  all  the  atoms  recorded  every 
picosecond.  The  first  100  ps  were  discarded  as  anneal  time.  Then  density  probability  functions 
p(z)  were  constructed  by  averaging  over  the  remaining  time  and  over  the  x  and  y  coordinates 
parallel  to  the  metal  surface  to  give  functions  of  z  the  coordinate  in  the  direction  perpendicular 
to  the  surface.  Figure  I,  in  particular  shows  the  density  profile  results  for  Na*^  cations.  Cl'  anions 
and  SPCE  water  molecules  in  a  3M  salt  solution  at  a  temperature  of  S0°  C.  Calculations  with 
temperatures  as  high  as  100°  C  where  originally  selected  to  improve  the  statistics.  Temperature 
dependence  of  the  probability  distributions  in  the  range  30  to  100°  C  is  weak.  In  the  50°  C 
simulation  there  were  94  Na  ions,  86  Cl  ions  and  1416  water  molecules.  The  charge  on  the 
electrode,  equal  to  the  image  charge,  is  -(94  -  86)e  or  -8e,  equivalent  to  a  surface  charge  density 
of -17.77  Cm-2. 

The  excess  sodium  ions,  attracted  to  the  electrode  by  the  image  charge,  screen  the  electrode 
charge.  The  screening  charge  is  concentrated  in  the  peak  at  1.45  nm  in  the  sodium  density  on 
the  right  hand  side  of  the  simulation  box  near  the  point  where  the  9-3  wall  potential  passes 
through  zero  (z^  =  1.615  nm).  The  image  plane  of  the  metal  is  located  at  the  extreme  the  right 
hand  vertical  edge  of  the  Figure  1  at  z  =  1.862  nm.  Simple  Gouy-(Zhapman  theory  would  predict 
an  increase  in  charge  to  a  finite  value  at  the  metal  surface.  There  is  no  peak  in  simple  theory 
because  ions  are  point-like  objects  and  the  solvent  is  a  continuum  with  a  fixed  dielectric  con¬ 
stant.  The  peak  in  the  cation  distribution  at  z  =  1.45  nm  is  a  consequence  of  a  fmite  ion  radius 
and  the  hydration  of  the  ions.  The  solvent  layer  at  the  electrode  also  effects  the  position  and 
shape  of  the  cation  distribution  near  the  metal.  Note  that  for  Izl  <  1.0  nm  the  cation  distribution 
is  approximately  flat  at  30  ions  nm'*.  There  is  also  a  small  peak  on  the  left  hand  side  at  ca.  z 
=  -1.4  nm,  that  is  not  associated  with  screening  but  is  instead  due  to  layering  of  the  water 
molecules  at  the  flat  restraining  wall.  This  wall  potential  goes  through  zero  at  z  =  -1.615  nm. 
Note  that  the  chloride  ion  distribution  does  not  show  a  peak  in  the  same  region  as  the  cation 
peak  because  negative  ions  are  repulsed  by  the  negatively  charged  metal  electrode. 

The  chloride  probability  distribution  has  no  major  structural  features,  certainly  no  peak  like  the 
screening  peak  at  z  =  1 .45  nm  in  the  Na''’  distribution.  Starting  from  the  metal  on  the  right  side 
of  Figure  1,  the  chloride  ion  distribution  rises  to  a  rough  plateau  that  stretches  from  z  =  1.00 
nm  to  z  =  -1.00  nm.  The  chloride  and  sodium  ion  probabilities  are  sufficiently  similar  across 
the  plateau  region  for  us  to  propose  that  this  has  the  properties  of  a  bulk  zone.  This  suggesticm 
is  supported  by  the  observation  of  approximate  local  charge  neutrality  (ion  densities  are  the 
same),  and  very  nearly  equal  integrated  densities  shown  in  Figure  2  for  each  of  the  concen¬ 
trations  at  which  the  calculations  were  performed.  The  integrated  density  plots  all  mcaiotcmically 
increase  with  z  from  left  to  right.  The  initial  slopes  for  z  <  -1.0  nm  are  less  because  of  the  rise 
In  ion  densities  from  zero  at  z  =  -1.615  nm  to  'bulk'  values  for  z  >  -1.0  nm.  The  solid  line  is 
the  integrated  density  for  the  Na'''  ion  and  the  dash  line  is  for  the  Q'  ions.  Since  the  integrated 
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densities  are  very  similar  this  means  that  the  electric  field  is  locally  very  small  or  zero  in  the 
bulk  zone  defined  by  Izl  <  1.0  nm. 

Whereas  on  the  basis  of  one  simulation  we  cannot  interpret  any  of  the  minor  features  in  the  ion 
profiles,  comparison  of  the  three  sets  in  Figure  2  shows  the  persistence  of  smaller  features  near 
1.2  nm  and  1.4  nm  for  both  ions.  We  will  discuss  diese  features  below. 

The  water  density  profile  shown  in  Figure  1  displays  well  pronounced  oscillations  on  both  sides 
of  the  simulation  cell  due  to  the  presence  of  the  walls  at  izl  =  1 .862  nm.  Recall  that  the  9-3 
wall  potential  is  the  same  on  both  sides.  The  metal  side  (right)  is  distinguished  by  the  existence 
of  the  image  plane.  On  the  metal  side  at  z  =  1.615  nm  the  water  density  peak  is  the  stnmgest 
feature  in  the  system  due  to  localized  water  next  to  the  metal.  The  surface  electric  field  reaches 
its  highest  value  in  this  region  and  water  is  localized  as  a  result  of  its  interacticMi  with  the  surface 
field.  In  this  simulation  the  waters  in  the  first  layer  are  oriented  by  this  strong  field  and  some 
H  bonds  to  the  bulk  region  are  broken.  This  electric  field  effect  is  distinct  from  localization  of 
water  cot  Pt(lOO)  and  Pt(lll)  surfaces  in  the  simulations  of  Heinzinger  and  Spohr'*,  and 
Berkowitz*^'  In  these  later  works  the  localization  is  at  top  sites  of  the  Pt  surface  due  to  di¬ 
rected  features  in  the  weak  chemisorptive  potential. 

In  Figure  2  we  show  the  results  of  calculations  that  show  the  dependence  of  the  ion  densities 
on  salt  concentration  and  electrode  charge.  The  30°  C  temperature  for  the  2M  and  IM  plots 
was  chosen  because  there  did  not  appear  to  be  a  significant  improvement  in  statistics  in  the 
earlier  50“  C  3M  calculation.  Each  calculation  took  approximately  3  months  on  a  work  station 
80  to  90%  dedicated  to  these  calculations  so  it  was  not  possible  to  run  all  the  calculations  at  the 
same  temperature.  In  these  plots  the  Cl  density  has  been  smoothed  to  distinguish  it  from  the 
Na  ion  density.  Top  panel  displays  the  3M  result  of  Figure  1  in  which  there  were  94  Na  ions, 
86  Cl  ions  and  1416  water  molecules  at  100“C,  and  the  image  charge  on  the  electrode  is  -<94  - 
86)e  =  -8e.  The  temperature  is  50“C.  Note  the  similarity  in  ion  densities  between  -1.0  and  1.0 
nm  in  all  three  examples.  The  middle  panel  displays  the  2M  result  in  which  there  were  62  Na 
ions,  58  Cl  ions  and  1516  water  molecules,  and  the  image  charge  on  the  electrode  is  -4e.  The 
temperature  is  30“C.  The  bottom  panel  displays  the  IM  result  in  which  there  were  32  Na  ions, 
28  Cl  ions  and  1576  water  molecules,  and  the  image  charge  on  the  electrode  is  again  -4e.  The 
temperature  is  30“  C. 

In  each  panel  the  cation  and  anion  concentrations  are  approximately  the  same  for  Izl  <1.0  nm. 
This  identifies  the  region  of  the  system  with  bulk  electrolyte  properties.  Also  for  each  system 
we  show  the  integrated  ion  density  vs  distance  z.  The  integrals  are  taken  from  z  =  -  1.862  nm 
to  z  =  1 .862  nm..  Solid  line  for  the  cation  and  dash  line  for  the  anion.  £)etailed  examination 
shows  that  in  each  case  for  Izl  <  1.0  nm  the  integrals  are  almost  indistinguishable.  This  indicates 
that  the  electric  field  due  to  the  ions  is  almost  zero  across  the  film  as  expected  for  a  bulk 
electrolyte. 

In  each  case  studied  the  width  of  the  electric  double  layer  region  is  small  but  finite.  As  an  es¬ 
timate  of  diffuse  layer  width  we  measure  from  the  position  of  the  hard  wall  (z  =  1.615  nm)  to 
the  point  where  the  difference  in  integrated  ion  densities  is  0.35  of  the  difference  at  the  metal 
(either  8  for  3M  or  4  for  IM  or  2M).  V/e  can  compare  this  with  the  Debye-Hiickel  screening 
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length  since  this  is  a  e''  measure  of  diffuse  layer  thickness  in  the  linearized  form  of  the  Gouy- 
Chapman  theory  (numbers  in  parentheses).  For  the  three  solutions  the  measures  are:  for  3M 
it  is  0.2  nm  (0.18  nm),  for  2M  it  is  0.4  nm  (0.24  nm),  for  I M  it  is  0.S  nm  (0.31  nm).  The 
numbers  in  parentheses  are  the  Debye  shielding  lengths  calculated  using  Eqn(l).  The  estimates 
based  on  the  ion  probability  distributions  must  be  regarded  as  very  rough  because  of  the  marmer 
they  were  done.  They  may  change  when  the  simulations  are  run  longer,  and  are  very  sensitive 
to  structure  in  the  water  layers  near  the  metal.  In  reality  there  must  be  effects  due  to  the  finite 
size  of  water  molecules  and  ions,  and  effects  due  to  saturation  of  molecular  polarizability  in  high 
electric  fields.  However  there  is  general  agreement  in  trend  and  magnitude.  This  agreement 
would  be  lost  if  the  dielectric  constant  for  water  were  taken  to  be  e  =  6,  a  value  more  appropriate 
for  a  zone  in  which  the  dielectric  properties  of  water  are  at  saturation  values*^. 

Returning  to  the  topic  of  fine  structure  in  the  ion  density  profiles  we  note  that  on  the  metal  side 
all  the  chloride  distributions  have  two  weak  features  at  ca.  1.2  and  1.4  nm.  Both  features  appear 
to  be  associated  with  peaks  in  the  cation  distribution  and  may  therefore  be  due  to  contact  pairs 
or  solvent  separated  pairs.  It  should  be  cautioned  that  the  SPC£  model  for  water  was  not  de¬ 
signed  with  high  salt  concentrations  in  mind,  so  the  ion  pairs  may  be  more  a  feature  of  the  model 
and  not  nature!  This  in  turn  implies  that  correlation  between  ions  at  high  salt  concentrations 
effects  the  distribution  near  the  charged  surface.  It  is  interesting  to  note  that  there  is  no  apparent 
association  of  similarly  positioned  ion  peaks  at  the  restraining  wall  on  the  left  hand  side  of 
Figure  2. 

Layering  in  the  water  distribution  effects  the  ion  distributions.  Figure  I  shows  a  general  result 
confirmed  by  more  extensive  simulations,  namely  ions  avoid  regions  of  high  water  density. 
On  the  metal  side  at  right  this  is  at  z  -  1.61S  nm  where  the  wall  potential  is  zero,  and  at  z  = 
1.25  nm  where  there  is  a  secondary  water  peak.  Note  that  the  highly  localized  water  layer  next 
to  the  electrode  creates  a  large  exclusion  region  for  ions  of  both  signs  compared  to  the  restraining 
wall  cn  the  left  hand  side  of  the  cell. 
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SUMMARY 


This  letter  described  the  calculation  of  die  structure  of  an  electric  double  layer  in  a  system  large 
enough  to  show  a  bulk-like  electrolyte  region  that  occupied  a  significant  portion  of  the  whole 
simulation  cell.  Estimates  of  the  width  of  the  diffuse  layer  were  consistently  larger  than  the 
corresponding  Debye-HUckel  screening  lengths. 

The  cations  and  anions  density  profiles  screening  the  electrode  was  structured.  In  the  case  of 
cations  diete  was  a  peak  near  1 .4  nm.  For  anions  there  was  a  shoulder  at  the  same  position. 
There  no  sign  of  similar  structure  near  the  restraining  wall.  Influence  of  solvent  on  both  dis¬ 
tributions  near  the  metal  interface  was  also  noted.  This  simulation  places  no  restrictions  on  ei¬ 
ther  ion-ion  or  water-ion  correlations,  unlike  the  simple  Gouy-Chapman  theory  where  there  is 
no  correlation. 

The  water  distribution  showed  a  plane  of  localized  water  at  the  point  where  the  9  -  3  wall  po¬ 
tential  passed  through  zero.  This  localization  involves  preferential  orientation  of  the  water,  and 
some  evidence  of  up  to  three  layers  of  water  is  evident. 
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Figure  1.  Ion  and  water  probability  density  distributions  for  3M  salt  solution  at  S0°C,  charge 
on  the  electrode  q,^  =  -8e.  Composition  of  the  simulation  cell:  94  Na*  ions,  86  Cr  ions,  and 
1463  water  molecules.  Metal  electrode  on  right  hand  side  at  z  =  1.862nm  (position  of  image 
plane),  restraining  wall  origin  on  the  left  at  z  =  -1.862  nm,  positions  where  the  wall  potentials 
on  both  sides  of  the  simulation  cell  go  through  zero  at  izi=  1.6IS  nm  are  shown  by  vertical 
broken  lines.  Water  density  shows  highly  localized  surf^  layer  of  water  near  z  =  1.62  nm. 
The  cation  distribution  shows  a  sharp  peak  near  z  »  1.4S  nm  corresponding  to  the  excess  positive 
ions  that  screen  the  negatively  charged  electrode.  The  Gouy-Chapman  distribution  for  cations 
and  anions  is  shown  superimposed  on  the  chloride  (broken  line)  and  sodium  (solid  line)  ion 
distribution. 
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Figure  2.  The  sodium  and  chloride  ion  probability  distributions  in  3M  (S0°C,  electrode  charge 
-8e),  2M  (30°C,  electrode  charge  -4e),  and  IM  (30®C,  electrode  charge  -4e)  aqueous  sodium 
chloride  solution.  Smoothed  curves  are  for  chloride,  and  the  unsmoothed  curves  are  for  sodium 
cations.  Metal  electrode  on  right  hand  side  at  z  =  !.862nm  (position  of  image  plane),  restraining 
wall  origin  on  the  left  at  z  =  -1.862  nm,  wall  potentials  on  both  sides  of  the  simulation  cell  go 
through  zero  at  lzl=  1.615  nm  (shown  by  vertical  broken  lines).  The  integral  of  each  ion  dis¬ 
tribution  is  plotted  ,  solid  line  for  cation  dashed  line  for  the  anion.  The  two  integrals  diverge 
near  the  metal  wall,  where  the  excess  cations  are  concentrated  and  shield  the  electrode  charge. 
The  depletion  of  the  anion  distribution  has  structure  that  is  induced  by  layering  in  the  water 
profile.  Note  that  all  the  distributions  show  a  distance  dependence  that  deviates  from  the  be¬ 
havior  expected  qualitatively  from  the  simple  Gouy-Chapman  theory. 
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